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ABSTRACT 

This paper presents a systematic treatment of the linear theory of scalar 
gravitational perturbations in the synchronous gauge and the conformal Newtonian 
(or longitudinal) gauge. It differs from others in the literature in that we give, in 
both gauges, a complete discussion of all particle species that are relevant to any flat 
cold dark matter (CDM), hot dark matter (HDM), or CDM+HDM models (including 
a possible cosmological constant). The particles considered include CDM, baryons, 
photons, massless neutrinos, and massive neutrinos (an HDM candidate), where the 
CDM and baryons are treated as fluids while a detailed phase-space description is 
given to the photons and neutrinos. Particular care is applied to the massive neutrino 
component, which has been either ignored or approximated crudely in previous works. 
Isentropic initial conditions on super-horizon scales are derived. The coupled, linearized 
Boltzmann, Einstein and fluid equations that govern the evolution of the metric and 
density perturbations are then solved numerically in both gauges for the standard 
CDM model and two CDM+HDM models with neutrino mass densities VL V = 0.2 and 
0.3, assuming a scale-invariant, adiabatic spectrum of primordial fluctuations. We 
also give the full details of the cosmic microwave background anisotropy, and present 
the first accurate calculations of the angular power spectra in the two CDM+HDM 
models including photon polarization, higher neutrino multipole moments, and 
helium recombination. The numerical programs for both gauges are available at 



http : //arcturus .mit . edu/cosmics 
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1. Introduction 

The theory of galaxy formation based on gravitational instability is aimed at describing how 
primordially-generated fluctuations in matter and radiation grow into galaxies and clusters of 
galaxies due to self-gravity. A perturbation theory can be formulated when the amplitudes of the 
fluctuations are small, and the growth of the fluctuations can be solved from the linear theory. 

In the early universe, gravitational perturbations are inflated to wavelengths beyond the 
horizon at the end of the inflationary epoch. Fluctuations of a given length scale reenter the 
horizon at a later time when the horizon has grown to the size of the fluctuations. Although 
the process of galaxy formation in recent epochs is well described by Newtonian gravity (and 
other physical processes such as hydrodynamics), a general relativistic treatment is required for 
perturbations on scales larger than the horizon size before the horizon crossing time. The use of 
general relativity brought in the issue of gauge freedom which has caused some confusion over 
the years. Lifshitz (1946) adopted the "synchronous gauge" for his coordinate system, which has 
since become the most commonly used gauge for cosmological perturbation theories. However, 
some complications associated with this gauge such as the appearance of coordinate singularities 
and spurious gauge modes prompted Bardeen (1980) and others (e.g. Kodama & Sasaki, 1984) 
to formulate alternative approaches that deal only with gauge-invariant quantities. A thorough 
review of the gauge-invariant perturbation theory and its application to texture-seeded structure 
formation models is given by Durrer (1993). Another possibility is to adopt a different gauge. 
We will discuss in detail in this paper the conformal Newtonian (or the longitudinal) gauge 
(Mukhanov, Feldman & Brandenberger 1992), which is a particularly convenient gauge to use for 
scalar perturbations. It also has the feature that the two scalar fields that describe the metric 
perturbations in this gauge are identical (up to a minus sign) to the gauge-invariant variables 
Bardeen (1980) constructed. 

The linear theory for a perturbed Friedmann-Robertson- Walker universe was first developed 
by Lifshitz (1946), later reviewed in Lifshitz & Khalatnikov (1963). The subsequent work can 
be found summarized in the textbooks by Weinberg (1972) and Peebles (1980), in the reviews 
by Kodama & Sasaki (1984) and Mukhanov et al (1992), and in the Summer School lectures 
by Efstathiou (1990), Bertschinger (1995), and Bond (1995). The theory has been applied to 
various cosmological models. In the synchronous gauge, for example, Peebles & Yu (1970) solved 
the linear evolution equations for photons and baryons in the absence of massless neutrinos and 
non-baryonic dark matter. Bond and Szalay (1983) included the neutrinos but approximated them 
as a perfect fluid by using only the first two moments in the angular expansion of the neutrino 
distribution function. They also did not include CDM. Xiang & Kiang (1992) considered both 
CDM and massive neutrinos, but treated the latter as non-relativistic particles and also ignored 
massless neutrinos. Holtzman (1989) provided fitting formulas for the baryonic transfer functions 
in various cosmological models, but gave no detailed description of the calculation, and again 
included only the first two moments of the neutrino distribution function. Durrer (1989) obtained 
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numerical solutions for the perturbations in the gauge-invariant formalism in the absence of 
baryons. Instead of using the convenient multipole expansion technique, Stompor (1994) adopted 
the earliest method by solving the massive neutrino distribution function directly at discrete 
momenta and angles. 

Another important application of the linear theory is to the cosmic microwave background 
anisotropies. The early calculations of the radiation perturbations included only the baryons 
and the photons (Peebles Sz Yu 1970; Wilson & Silk 1981). The calculations were subsequently 
extended to dark matter dominated models (e.g. Bond & Efstathiou 1984, 1987; Vittorio & Silk 
1984, 1992; Holtzman 1989). The post-COBE era has seen a rapid increase in the number of 
papers on this subject (see Sugiyama & Gouda 1992; White, Scott, & Silk 1994; and references 
therein). However, there has not been a complete treatment of the massive neutrinos to the 
accuracy that is needed for comparison with the near-future anisotropy experiments. Truncating 
the massless neutrino moments beyond I > 2 in the pure CDM model leads to a 10% error in the 
anisotropy power spectrum (Hu et al 1995). One would expect a comparable or larger error if the 
I > 2 modes were not included for the massless and massive neutrinos in CDM+HDM models. 

This paper differs from earlier ones in the literature in that we give a complete discussion of 
all particle species that are relevant to any flat CDM, HDM, or CDM+HDM models (including 
a possible cosmological constant). These include CDM, baryons, photons, massless neutrinos, 
and massive neutrinos (an HDM candidate). The recent interests in CDM+HDM models (e.g. 
Davis, Summers, & Schlegel 1992; Klypin et al 1993; Jing et al 1994; Cen k Ostriker 1994; Ma 
& Bertschinger 1994b; Ma 1994; Primack et al 1995; Klypin et al 1995) prompted us to provide 
a detailed treatment of the massive neutrino component, which was either ignored or crudely 
approximated in other works. Part of the cause for neglecting this component was perhaps a lack 
of motivation when the standard CDM model was thought to match the observed structure well. 
Another was perhaps due the sharp increase in the computing time required to deal with the 
time-dependent nature of the massive neutrinos' energy-momentum relation. For completeness, we 
also give full details and results of the calculation of the cosmic microwave background anisotropy 
including photon polarization and helium recombination. 

This paper serves two purposes. First, it is an independent paper in which we present a 
systematic treatment of the linear theory of scalar isentropic gravitational perturbations in the 
synchronous and conformal Newtonian gauges. The coupled, linearized Einstein, Boltzmann, and 
fluid equations for the metric and density perturbations are presented in parallel in the two gauges. 
The CDM and the baryon components behave like collisionless and collisional fluids, respectively, 
while the photons and the neutrinos require a phase-space description governed by the Boltzmann 
transport equation. We also derive analytically the time dependence of the perturbations on 
scales larger than the horizon to illustrate the dependence on the gauge choice. This information 
is needed in the initial conditions for the numerical integration of the evolution equations. 

This paper also serves as a companion paper to Ma & Bertschinger (1994a) in which we 
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reported the main results from our linear calculation of the full neutrino phase space in a 
CDM+HDM model with ft c = 0.65, Q v = 0.3, O baryon = 0.05, and H = 50 km s" 1 Mpc" 1 . (The 
corresponding neutrino mass is m u = 7 eV.) The motivation was to obtain an accurate sampling 
of the neutrino phase space for the HDM initial conditions in subsequent iV-body simulations 
(Ma & Bertschinger 1994b; Ma 1994) of structure formation in CDM+HDM models. We adopted 
a two-step Monte Carlo procedure to achieve this goal: (1) Integrate the coupled, linearized 
Boltzmann, Einstein, and fluid equations for all particle species in the model (i.e., CDM, HDM, 
photons, baryons and massless neutrinos) to obtain the evolution of the metric perturbations; (2) 
Follow the trajectories of individual neutrinos by integrating the geodesic equations using the 
metric computed in (1). Since no coordinate singularities occur in the conformal Newtonian gauge 
and the geodesic equations have simple forms, the geodesic integration in step (2) was carried out 
in this gauge, starting shortly after neutrino decoupling at redshift z ~ 10 9 until z = 13.5. We 
focus on step (1) in this paper. Following historical precedents, we first developed the code for the 
Boltzmann integration in the synchronous gauge. The transformation relating the synchronous 
gauge and the conformal Newtonian gauge was then derived and used to compute the metric 
perturbations in the latter gauge for step (2) of the calculation. Subsequently we developed a code 
to perform the full integration in the conformal Newtonian gauge. 

The organization of this paper is as follows. In §2 we write down the metric for the two gauges 
and summarize their properties. In §3, we derive the gauge transformation relating two arbitrary 
gauges and obtain the transformation between the synchronous and the conformal Newtonian 
gauges. The linearized evolution equations for the metric and the density perturbations are given 
in §§4 and 5. Section 4 discusses the Einstein equations with emphasis on the source terms, the 
energy-momentum tensor, in the two gauges. The perturbed fluid equations are derived from the 
energy-momentum conservation, which are applied to CDM and the baryons in §5. The rest of §5 
contains detailed treatments of the photon and neutrino phase space distributions, recombination, 
and the coupling of photons and baryons. The photon and neutrino distribution functions are 
expanded in Legendre polynomials, reducing the linearized Boltzmann equation to a set of coupled 
ordinary differential equations for the expansion modes. The massive neutrinos require a slightly 
more complicated treatment because the energy-momentum relation depends nontrivially on time 
when the neutrinos make a transition from being relativistic to non-relativistic. Our method 
for computing the microwave background anisotropy is presented in §6. Section 7 discusses the 
behavior of the perturbations before horizon crossing. The necessary initial conditions for the 
variables in the two gauges are given. Section 8 presents the numerical results for the evolution of 
the perturbations and the angular power spectrum of the microwave anisotropy in CDM+HDM 
models. 

The complex physics we study requires the use of many equations and symbols. As a guide 
to the reader, in Table 1 we summarize the key symbols and the equations where they are defined 
or first used. 
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Symbol 


Meaning 


Reference Eqn 


a 


scale factor 


(1) 


T 


conformal time 


(1) 


k 


wavenumber of Fourier mode 


© 


Pi 


conjugate momentum to (comoving) position x l 


(31) 


Pi 


proper momentum 


dD 


qi 


= aPi 




e 


= (q 2 + a 2 m 2 ) 1 ' 2 


— 


h 


synchronous 


a), & 


V 


metric perturbations 






conformal Newtonian 


© 




metric perturbations 




subscript c 


cold dark matter 


- 


subscript v 


massless neutrinos 


— 


subscript h 


massive neutrinos 


- 


subscript 7 


photons 


- 


subscript b 


baryons 


- 


fo 


unperturbed phase space distribution function 


dD 




perturbation to fo 






1th Legendre component of 


© 


Fi 


momentum-averaged 


dD 




photon polarization component 


dH) 


S 


density fluctuation (= Fq) 


(Ha), © 


e 


divergence of fluid velocity (= 3/cFi/4) 


(22), @>), © 


a 


shear stress (= F2/2) 


(H), &), © 


A 


temperature fluctuation (= AT/T = F 7 /4) 


dD 


A/ 


Zth Legendre mode of A (= F 7 //4) 


(84) 


Q 


temperature fluctuation power spectrum 


© 


C s 


baryon sound speed (= SP^/Sp^) 




Cpb 


sound speed of coupled photon-baryon fluid 


(H) 


w 


describes equation of state (= P/p) 




n e 


electron number density 


© 




Thomson scattering cross section 


dD 


T c 


= (an e CTT) -1 


dD 


Ry 


density ratio pv/ip^ + p v ) 


dD 


R 


density ratio (4/3)p 7 /ph 


diO) 



Table 1: Symbols used in this paper. 
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2. The Two Gauges 

We consider only spatially flat (Q = 1) background spacetimes with isentropic scalar metric 
perturbations. The spacetime coordinates are denoted by x M , n G (0,1,2,3), where x° is the 
time component and x l , i £ (1,2,3) are the spatial components in Cartesian coordinates. Greek 
letters a, (3, 7 and so on always run from to 3, labeling the four spacetime-coordinates; Roman 
letters such as i,j, k always run from 1 to 3, labeling the spatial parts of a four-vector. Repeated 
indices are summed. Since our interests lie in the physics in an expanding universe, we use 
comoving coordinates x M = (r, x) with the expansion factor a(r) of the universe factored out. The 
comoving coordinates are related to the proper time and positions t and r by dx° = dr = dt/a(r), 
dx = dr/a(r). Dots will denote derivatives with respect to r: a = da/dr. The speed of light c is 
set to unity. 

The components goo and go% of the metric tensor in the synchronous gauge are by definition 
unperturbed. The line element is given by 

ds 2 = a 2 {r){-dT 2 + (Sij + hifidjdx?} . (1) 

The metric perturbation hij can be decomposed into a trace part h = ha and a traceless part 
consisting of three pieces, hfj, hjj, and hfj, where = hSij/3 + hjj + hjj + hfj. By definition, the 
divergences of h}\- and hjj (which are vectors) are longitudinal and transverse, respectively, and hfj 
is transverse, satisfying 

eijkdjdih)} k = , didjh- - = , dihfj = . (2) 
It then follows that h)) H can be written in terms of some scalar field u and hj~, in terms of some 

l 3 ' L J 

divergenceless vector A as 

hij = (didj ~ ^<% v2 ) ^ ' 

h± = fyAj+djAi, d t A i = 0. (3) 

The two scalar fields h and \i (or h\-) characterize the scalar mode of the metric perturbations, 
while Ai (or hjj) and hfj represent the vector and the tensor modes, respectively. 

We will be working in the Fourier space k in this paper. We introduce two fields h(k, r) and 
rj(k,r) in /c-space and write the scalar mode of hij as a Fourier integral 

hij(x, t) = J d 3 ke^' s ^kikjh(k, r) + (kfkj - ^^) 6rj(k, r)| , k = kk . (4) 

Note that h is used to denote the trace of in both the real space and the Fourier space. 

In spite of its wide-spread use, there are serious disadvantages associated with the synchronous 
gauge. Since the choice of the initial hypersurface and its coordinate assignments are arbitrary, 
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the synchronous gauge conditions do not fix the gauge degrees of freedom completely. Such 
residual gauge freedom is manifested in the spurious gauge modes contained in the solutions to 
the equations for the density perturbations. The appearance of these modes has caused some 
confusion over the years and prompted Bardeen (1980) to formulate alternative approaches that 
deal only with gauge- invariant quantities. Another difficulty with the synchronous gauge is that 
since the coordinates are defined by freely falling observers, coordinate singularities arise when two 
observers' trajectories intersect each other: a point in spacetime will have two coordinate labels. 
A different initial hypersurface of constant time has to be chosen to remove these singularities. 

The conformal Newtonian gauge (also known as the longitudinal gauge) advocated by 
Mukhanov et al. (1992) is a particularly simple gauge to use for the scalar mode of metric 
perturbations. The perturbations are characterized by two scalar potentials ip and (ft which appear 
in the line element as 

ds 2 = a 2 (r) {-(1 + 2^)dr 2 + (1 - 2cft)dx i dx i \ . (5) 

It should be emphasized that the conformal Newtonian gauge is a restricted gauge since the metric 
is applicable only for the scalar mode of the metric perturbations; the vector and the tensor 
degrees of freedom are eliminated from the beginning. Nonetheless, it can be easily generalized to 
include the vector and the tensor modes (Bertschinger 1995). The discussion here will be confined 
to the scalar perturbations only. 

One advantage of working in this gauge is that the metric tensor g^ v is diagonal. This 
simplifies the calculations and leads to simple geodesic equations (Ma &: Bertschinger 1994a). 
Another advantage is that ip plays the role of the gravitational potential in the Newtonian limit 
and thus has a simple physical interpretation. Moreover, the two scalar potentials ip and (ft in this 
gauge are related to the gauge-invariant variables §a and of Bardeen (1980) and \E' and $ of 
Kodama & Sasaki (1984) by the simple relation 

iP = $ A = ifr } (ft=-$ H = -$. (6) 

No gauge modes are present in this gauge to obscure the meaning of the physical modes since 
the gauge freedom is completely fixed for 0, = 1 aside from the addition of spatial constants to ift 
and (ft. The two potentials differ when the energy-momentum tensor T^ v contains a nonvanishing 
traceless and longitudinal component (see eq. [|23d|1 ). 



3. Gauge Transformations 

In this section we first derive the transformation law relating two arbitrary gauges. From it, 
the gauge transformation relating the synchronous gauge and the conformal Newtonian gauge is 
readily obtained. 
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A perturbed flat Friedmann-Robertson- Walker metric can be written in general as 
5oo = -a 2 (r) {l + 2V(x,r)} , 

g 0i = a 2 (T)wi(x,T) , (7) 
gij = a 2 (r) {[1 - 24>(x,T)]5ij + Xij(x,r)} , xu = 

where the functions ip, </>, Wi and Xij represent metric perturbations about the Robertson- Walker 
spacetime and are assumed to be small compared with unity. The trace part of the perturbation 
to gij is absorbed in cfi, and Xij is taken to be traceless. 

Consider a general coordinate transformation from a coordinate system to another x^ 

x » _> x» = x» + d»(x u ) . (8) 
We write the time and the spatial parts separately as 

x° = x° + a(x, t) , 

f = x + V/3(f,r) +e(f,r), V-e = 0, (9) 

where the vector d has been decomposed into a longitudinal component V/3 (V x V/3 = 0) and a 
transverse component e (V • e = 0). The requirement that ds 2 be invariant under this coordinate 
transformation leads to 

9nu{x) = 9ixu(x) - g fl p(x)d u d 13 - g a v{x)d PL d a - d a d a g fiu (x) + 0(d 2 ) . (10) 

We note that both sides of this equation are evaluated at the same coordinate values x in the two 
gauges, which do not correspond to the same physical point in general. Assuming to be of 
the same order as the metric perturbations ip,Wi,4> an d Xiji the metric perturbations in the two 
coordinate systems are related to first order in the perturbed quantities by 

^(x,t) = t/j(x,t) — a(x,r) a(x, r) , (H a ) 

Wi(x,T) = Wi(x,T)+diCt(x,T) -diP(x,T)-€i(x,T) , (lib) 
4>(x,t) = (j)(x,T) + \v 2 l3{x,T) + -a(x,T), (11c) 

O Ob 

Xij(x,r) = Xi ,0r,T)-2{(^-^%^ (lid) 

We can further decompose the transformations of Wi and Xij above into longitudinal and 
transverse parts: 

w\{x,t) = wj(x,T) + dia(x,r) - di$(x,T) , 

wI{x,t) = wI{x,t) - ei{x,r) , (12) 



- 9 - 



and 



2 ( d l d ] - -5 tl V 2 



/3(x, t) 



Xij (x, t) = xh (x, t) - {diEj + djei) , 



Xij i x i 7 " 



Xij ( x ? T ) 



(13) 



where Wi 



w i + Wi, Xij = Xij + Xij + Xij , and x\j , Xij and xfj obey equation (|). Equations 



(11 )— (13) describe the transformation of metric perturbations under a general infinitesimal 
coordinate transformation. 



We can now use equations ( |ll| ) to relate the scalar metric perturbations (0, ip) in the 
conformal Newtonian gauge to hij = h5ij/3 + h)\- in the synchronous gauge. Let x^ 1 denote the 
synchronous coordinates and the conformal Newtonian coordinates with x^ = x^ + cP. From 
equations ( |i"2| ) and dl3|), we find 



a(x, r) 
ej(x,r) 

hjj(x,T) 

diCj + d,e 



/?(x,r)+£(r), 
e*(x) , 

1 

3"v 



-2 [ $9,- 



<5,;,V 2 



/3(x,r) 



0, 



(14a) 
(14b) 

(14c) 

(14d) 



where £(t) is an arbitrary function of time, reflecting the gauge freedom associated with the 
coordinate transformation: x° = x° + £(t), x l = x l . This transformation corresponds to a global 
redefinition of the units of time with no physical significance; therefore we shall set £ = from 
now on. From equations ( |lla| ) and ( |llc| ) we then obtain 



(15) 



V>(x,t) 
4>(x,t) 

and (3 is determined by /i" in equation ( |14cj ). 
In terms of h and rj introduced in equation 



+P(x, t) + -P(x , t) , 
a 



h(x,T)-\v 2 f3(x,T)--f3(x,r) 
bo a 



h-j in the synchronous gauge is given by 



hjj(x , t) = J d 3 k j ts {kikj - ^<%) [h{k, t) + 6r](k, r)} , k = kk . 



(16) 



Comparing h}\ - in equations (14c) and (16), we can read off (3: 



J ^ ke ^W {Hk,r) + Qn(k,r)} . 



(17) 



Then from equations ([15]), the conformal Newtonian potentials 4> and ifi are related to the 
synchronous potentials h and r\ in fc-space by 

^(k> r ) = 7^(/i(A?,r)+6f?(A?,r) + ^ [/i(fc,r)+67)(A?,r) 



2/fc 2 
rj(k, t) 



' r /i(fe,r) + 6rj(k,r) 



2k 2 



(18) 
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The other components of the metric perturbations, u>i, xtj , and xjj ; are zero i n both gauges. 



4. Einstein Equations and Energy-Momentum Conservation 

For a homogeneous Friedmann-Robertson- Walker universe with energy density p(r) and 
pressure P(r), the Einstein equations give the following evolution equations for the expansion 
factor a(r): 

-Y = %G*p-K, (19) 
a J 3 

= -%Ga\p + ZP), (20) 

where the dots denote derivatives with respect to r, and k is positive, zero, or negative for a 
closed, flat, or open universe, respectively. We consider only models with total 0, = 1 in this 
paper, so we set k = 0. A cosmological constant is allowed through its inclusion in p and P: 
PA = A/87rG = — Pa- This is the only place that A enters in the entire set of calculations. It follows 
from equation (^) with k = that the expansion factor scales as a oc r in the radiation-dominated 
era, a oc r 2 in the matter-dominated era, and a oc (too — r ) _1 m a cosmological constant-dominated 
era (in the latter the radius of the de Sitter event horizon). 

We find it most convenient to solve the linearized Einstein equations in the two gauges in 
the Fourier space k. In the synchronous gauge, the scalar perturbations are characterized by 
h(k,r) and rf(k,r) in equation In terms of h and rj, the time-time, longitudinal time-space, 
trace space-space, and longitudinal traceless space-space parts of the Einstein equations give the 
following four equations to linear order in A;-space: 

Synchronous gauge — 

k 2 r,---h = 4vrGa 2 5T o(Syn) , (21a) 
2 a 

k 2 r\ = 4vrGa 2 (p + P)^(Syn), (21b) 

h + 2-h-2k 2 rj = -8vrGa 2 <5TUSyn) , (21c) 
a 

h + 6rj + 2- f h + 67)) - 2k 2 n = -24vrGa 2 (p + P)a(Syn) . (21d) 
a V ' 

The label "Syn" is used to distinguish the components of the energy-momentum tensor in the 
synchronous gauge from those in the conformal Newtonian gauge. The variables and a are 
defined as 

(p + p)0 = ikidT ! , (p + P)a = -(kk - -Sijpj , (22) 

and X l j = T l j — S^T \/3 denotes the traceless component of T Kodama &: Sasaki (1984) define 
the anisotropic stress perturbation II, related to our a by a = 2IIP/3(p + P). When the different 
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components of matter and radiation (i.e., CDM, HDM, baryons, photons, and massless neutrinos) 
are treated separately, [p + P)6 = J2i(f>i + P~i)@i an d (p + P)o~ = J2i(f>i + Pi)°~i > where the index i 
runs over the particle species. 

In the conformal Newtonian gauge, the first-order perturbed Einstein equations give 
Conformal Newtonian gauge — 

a 



k 2 (p + 3 : 



k 





a J 


= 4vrGa 2 5T°o(Con) , 


(23a) 




U) 

a J 


= 4vrGa 2 (p + J P)0(Con), 


(23b) 


k\ 

— (6 
3 yv 




= —Ga 2 5T\(Con) , 


(23c) 


k\<P 


-1>) 


= 12TrGa 2 (p + P)a(Con), 


(23d) 



a , • •,. (a a? 

-(^ + 20+ 2 2 

a \ a a z 



where "Con" labels the conformal Newtonian coordinates. 

Now we derive the transformation relating 5T^ U in the two gauges. For a perfect fluid of 
energy density p and pressure P, the energy-momentum tensor has the form 

T% = Pg% + { P + P)U»U V , (24) 

where = dx^ /V —ds 2 is the four- velocity of the fluid. The pressure P and energy density p of 
a perfect fluid at a given point are defined to be the pressure and energy density measured by a 
comoving observer at rest with the fluid at the instant of measurements. For a fluid moving with 
a small coordinate velocity v % = dx l /dr, v % can be treated as a perturbation of the same order as 
bp = p — p, 5P = P — P, and the metric perturbations. Then to linear order in the perturbations 
the energy-momentum tensor is given by 

T° = -(p + Sp), 

T°i = (p + P) Vi = -T i , 

Tj = (P + 8P)5 i j + S l j , S\ = 0, (25) 

where we have allowed an anisotropic shear perturbation Ej in T l j. As we shall see, since the 
photons are tightly coupled to the baryons before recombination, the dominant contribution to 
this shear stress comes from the neutrinos. We note that for a fluid, 6 defined in equation (|22| ) is 
simply the divergence of the fluid velocity: 9 = ik 3 Vj. 

The energy-momentum tensor T'V(Syn) in the synchronous gauge is related to T^(Con) in 
the conformal Newtonian gauge by the transformation 

T ^ (Syn) = ^5F T ^ (COn) ' (26) 

where and denote the synchronous and the conformal Newtonian coordinates, respectively. 
It follows that to linear order, T° (Syn) = T° (Con) , T°j(Syn) = T°j(Con) + ikja(p + P) , and 
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T j(Syn) = T j(Con) , where a = x° — x° = (h + 6r])/2k 2 in fc-space from equations ( [L4a| ) and (|l7|). 
Let 5 = Sp/p = —5T%/p. Evaluating the perturbations at the same spacetime coordinate values, 
we obtain 



<5(Syn) = <5(Con)-c^, (27a) 

#(Syn) = ^(Con)-aA; 2 , (27b) 

<5P(Syn) = <5P(Con) - aP , (27c) 

a(Syn) = cr(Con) . (27d) 



This transformation also applies to individual species when more than one particle species 
contributes to the energy-momentum tensor, provided that the appropriate p and P are used for 
each component. 

The non-relativistic fluid description is appropriate for the CDM and the baryon components. 
The photon and the neutrino components, however, can be appropriately described only by their 
full distribution functions in phase space. The energy- momentum tensor in this case is expressed 
through integrals over momenta of the distribution functions. We will discuss it in detail in §5. 

The conservation of energy-momentum is a consequence of the Einstein equations. Let 
w = P/p describe the equation of state. Then the perturbed part of energy-momentum 
conservation equations 

= + T v ap T af} + T a ap T u P = (28) 

in A;-space implies 
Synchronous gauge — 

e = -°(i-3 ro )9-^e + ^*^i-iV, (29) 

a 1 + W 1 + W 



Conformal Newtonian gauge — 

= -^(l-3 w )e--^9+ S -^k 2 5-k 2 a + k 2 ^. (30) 
a 1+w 1+w 

These equations are valid for a single uncoupled fluid, or for the net (mass-averaged) 5 and 9 
for all fluids. They need to be modified for individual components if the components interact 
with each other. An example is the baryonic fluid in our model, which couples to the photons 
before recombination via Thomson scattering. In the next section we will show that an extra term 
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representing momentum transfer between the two components needs to be added to the 5 equation 
for the baryons. 

For the isentropic primordial perturbations considered in this paper, the equations above 
simplify since 5P = c 2 5p, where c 2 s = dP/dp = w + pdw/dp is the adiabatic sound speed squared. 
For the photons and baryons (the only collisional fluid components with pressure), w is a constant 
(w = 1/3 for photons and w ~ for baryons since they are nonrelativistic at the times of interest). 
Thus, 5P/5p — w = 0. The entropy generated from the coupling of photons and baryons before 
recombination is accounted for by Thomson scattering terms added to their respective equations 
of motion in §§5.5-5.7. Even in the case of isocurvature baryon models, which may have large 
perturbations in the entropy per baryon ab initio, the corrections to 5P/5p = w are generally very 
small because of the large photon to baryon ratio. 



5. Evolution Equations for Matter and Radiation 

5.1. Phase Space and the Boltzmann Equation 

A phase space is described by six variables: three positions x % and their conjugate momenta 
Pi. Our treatment of phase space is based on the time-slicing of a definite gauge (synchronous 
or conformal Newtonian). Although this approach is not manifestly covariant, it yields correct 
results provided the gauge-dependent quantities are converted to observables at the end of the 
computation. 

The conjugate momentum has the property that it is simply the spatial part of the 
4-momentum with lower indices, i.e., for a particle of mass m, Pi = mUi , where C7j = dxij \J —ds 1 . 
One can verify that the conjugate momentum is related to the proper momentum p % = pi measured 
by an observer at a fixed spatial coordinate value by 

Pi = a(5ij + ^hij)p> , in synchronous gauge , 

Pi = a(l — 4>)pi , in conformal Newtonian gauge . (31) 

In the absence of metric perturbations, Hamilton's equations imply that the conjugate momenta 
are constant, so the proper momenta redshift as a -1 . 

The phase space distribution of the particles gives the number of particles in a differential 
volume dx 1 dx 2 dx 3 dPidP 2 dP 3 in phase space: 

f(x\ P j ,T)dx 1 dx 2 dx 3 dPidP 2 dP 3 = dN . (32) 

Importantly, / is a scalar and is invariant under canonical transformations. The zeroth-order phase 
space distribution is the Fermi-Dirac distribution for fermions (+ sign) and the Bose-Einstein 
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distribution for bosons (— sign): 



/o = /o(e) = ^ e£/feBTo±1 , (33) 



where e = a(p 2 + to 2 ) 1 / 2 = (P 2 + a?m 2 ) 1 / 2 , To = aT denotes the temperature of the particles 
today, the factor g s is the number of spin degrees of freedom, and hp and /cb are the Planck and 
the Boltzmann constants. 

When the spacetime is perturbed, x % and P% remain canonically conjugate variables, with 
equations of motion given by Hamilton's equations (Bertschinger 1993). However, following 
common practice (e.g., Bond & Szalay 1983) we shall find it convenient to replace Pj by qj = apj 
in order to eliminate the metric perturbations from the definition of the momenta. Moreover, 
we shall write the comoving 3-momentum qj in terms of its magnitude and direction: qj = qnj 
where n l rii = bij-rtw? = 1. Thus, we change our phase space variables, replacing f(x l ,Pj,r) by 
f(x l , q, nj,r). While this is not a canonical transformation (i.e., qi is not the momentum conjugate 
to x' 1 ), it is perfectly valid provided that we correctly transform the momenta in Hamilton's 
equations. Note that we do not transform /. Because qj are not the conjugate momenta, d 3 xd 3 q 
is not the phase space volume element, and fd 3 xd 3 q is not the particle number. In the conformal 
Newtonian gauge, for example, (1 — 3(/))fd 3 xd 3 q is the particle number; this result is sensible 
because a(l — 4>)dx l is the proper distance. 

In the perturbed case we shall continue to define e as a(r) times the proper energy measured 
by a comoving observer, e = (q 2 + a 2 m 2 ) 1 / 2 . This is related to the time component of the 
4-momentum by Po = ~ e m the synchronous gauge and Po = — (1 + V0 e m the conformal 
Newtonian gauge. For the CDM+HDM models we are interested in, the photons, the massless 
neutrinos, and the massive neutrinos at the time of neutrino decoupling are all ultra-relativistic 
particles, so e in the unperturbed Fermi-Dirac and Bose-Einstein distributions can be simply 
replaced by the new variable q. 

The general expression for the energy-momentum tensor written in terms of the distribution 
function and the 4-momentum components is given by 

T, u = J dPxdP 2 dP 3 {-g)- 1 ' 2 ^f( x i, Pj , T ) , (34) 

where g denotes the determinant of g^ v . It is convenient to write the phase space distribution as a 
zeroth-order distribution plus a perturbed piece in the new variables q and ny. 

f(x\P j ,T)=Mq)\l + ^(x\q,n j ,r)} . (35) 



In the synchronous gauge, we have (— g)~ 1 ^ 2 = a~ 4 (l — ^h) and dPidP2dP% = (1 + ^h)q 2 dqdQ 
to linear order, where h = ha and d£l is the solid angle associated with direction rij. Using the 
relations / dQniUj = 4ir5ij/3 and / dVLrii = J dflninjiik = 0, it then follows from equation (|3~I| ) 
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that 

T°o = -a" 4 J q 2 dqdn \j q 2 + m 2 a 2 f (q) (1 + \&) , 

T°i = a" 4 J q 2 dqdnqn i f (q)^ , (36) 
J y'q z + m z a z 

to linear order in the perturbations. Note that we have eliminated the explicit dependence on the 
metric perturbations in equation (|34| ) by redefining Pi in terms of q and n{. Note also that the 
comoving energy e(q,r) = (q 2 + a 2 ™ 2 ) 1 / 2 is used in the integrands but not in the argument of the 
unperturbed distribution function. 

In the conformal Newtonian gauge, (— g)~ 1 ^ 2 = a _4 (l — ifi + 30) and dP\dP2dP^ = 
(1 — 3(j))q 2 dqd£l. It then follows that the components of the energy- momentum tensor have the 
same form as in equations (|^). Of course, it is understood that the variables q and rtj in this case 
are defined in relation to the conjugate momentum Pi in the conformal Newtonian coordinates, and 
not the synchronous coordinates. (They differ because comoving observers in the two coordinate 
systems are not the same.) The expansion factor a and ^ are evaluated at the coordinates (x l ,r) 
in the conformal Newtonian gauge. 

The phase space distribution evolves according to the Boltzmann equation. In terms of our 
variables (x' l ,q,nj,r) this is 

I^ = ^ + ^dJ_ + dqdl + dn 1 dl = /dl\ 
dr dr dr dx l dr dq dr dm \dr J q ' 

where the right-hand side involves terms due to collisions, whose form depends on the type of 
particle interactions involved. From the geodesic equation 



dP^ 

P°^- + P a P P = , (38) 



it is straightforward to show that 



dq/dr = — -qhijUiUj in synchronous gauge , 

dq/dr = q<p — e(q,r) mdiip in conformal Newtonian gauge , (39) 

and drii/dr is 0(h). Since df/drii is also a first-order quantity, the term (drii/dr)(df /dm) in the 
Boltzmann equation can be neglected to first order. Then the Boltzmann equation in fc-space can 
be written as 



Synchronous gauge 



— — h i -(k ■ n)* + — 

or e dmq 



^L(k.n) 2 



1 









(40) 

c 
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Conformal Newtonian gauge — 

rfinf„ r. e ^ -\ 1 /ftf \ 

(41) 



— h % - (k ■ n) V + 



dr e dlnq 



i —{% ■ n) ip 



1 









c 



Equations @ and @ can also be derived using the canonical phase space variables a;* and Pj 
and Hamilton's equations (instead of the geodesic equation), followed by a transformation from Pj 
to qrij. 

The terms in the Boltzmann equation depend on the direction of the momentum n only 
through its angle with k. (We shall see that this is true of the collision term for photons as well 
as the convective and metric perturbation terms.) Therefore, if the momentum-dependence of the 
initial phase space perturbation is axially symmetric about k, it will remain axially symmetric. 
If axially-asymmetric perturbations in the neutrinos or other collisionless particles are produced, 
they would generate no scalar metric perturbations and thus would have no effect on other species. 
Therefore, we shall assume that the initial momentum-dependence is axially symmetric so that 
depends on q = qh only through q and k ■ h. This assumption, which effectively reduces the 
dimensionality of phase space perturbations by one (after Fourier transforming on the spatial 
coordinates), has been made (implicitly, if not explicitly) in all previous studies of the evolution of 
scalar perturbations. 



5.2. Cold Dark Matter 



CDM interacts with other particles only through gravity and can be treated as a pressureless 
perfect fluid. The CDM particles can be used to define the synchronous coordinates and therefore 
have zero peculiar velocities in this gauge. Setting 8 = a = and w = w = in equations ( p9|) 
leads to 

Synchronous gauge — 

8c = -~h. (42) 

The CDM fluid velocity in the conformal Newtonian gauge, however, is not zero in general. In 
fc-space, equations (|30|) give 

Conformal Newtonian gauge — 

S c = -e c + 36, 9 C = --e c + k 2 ifj. (43) 

a 

The subscript c in 5 C and 9 C denotes the cold dark matter. 
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5.3. Massless Neutrinos 



The energy density and the pressure for massless neutrinos (labeled by subscripts v) are 
p v = 3P V = — T°o = T\. From equations (|36|) the unperturbed energy density p v and pressure P v 
are given by 

p v = 3P U = a- A J q 2 dqdn qf (q) , (44) 

and the perturbations of energy density 5p u , pressure 5P U7 energy flux ST^ , and shear stress 
T, l uj = T % vj - P u 5ij are given by 

5p u = 35P„ = a- A J q 2 dqdnqfo(q)*, 

5T° ui = a" 4 J q 2 dqdnqnif Q (q)^, (45) 

T, l vj = a~ 4 J q 2 dqdn q^Uj - -Sij) fo(q)V . 

The unperturbed energy flux and shear stress are zero. 

The Boltzmann equation simplifies for massless particles, for which e = q. To reduce the 
number of variables we integrate out the (/-dependence in the neutrino distribution function and 
expand the angular dependence of the perturbation in a series of Legendre polynomials Pi(k ■ n): 

F v (k, h, t) = J f^ qf °}f* = £H)'(2Z + l)F ul (k, r)P z (£ • h) . (46) 
Jq 2 dqqfo{q) ^ 

As noted in §5.1, the dependence on h arises only through k ■ h, so that a general distribution may 
be represented as in equation (f46[). The factor (— i) l (2l + 1) is chosen to simplify the expansion of 
a plane wave: F v = exp(—ik ■ x) with x = r(r)n has expansion coefficients F v \ = ji(kr) given by 
spherical Bessel functions. 

In terms of the new variable F u (k,h,T) and its harmonic expansion coefficients, the 
perturbations 5 V , p v , Q u , and a v (defined in eq. [p22| ) take the form 



- j dnF l/ (k,n,T)=F u0 , 

J dn(k-h)F u (k,n,T) = ^kF vl , (47) 



16tt 
3 

a v = / dU 

16vr 



{k-h) 2 - 1 - 



F u (k,h,r) = -F v2 



where we have used p u = 3P U for the massless neutrinos. 

Integrating equations (f40|) and (|4l| ) over q 2 dq qfo(q) and dividing them by / q 2 dqqfo(q), the 
Boltzmann equation for massless neutrinos becomes 
QP 2 ■ 4 ■ 

— h ikpF v = ~-h — — (h + Qrf)P2{p) in synchronous gauge , 

9F U 

— h ikpF v = 4 {4> — ikpvp) in conformal Newtonian gauge , (48) 



-18- 



where \i = k-n and Pi{p) = ^(3/i 2 — 1) is the Legendre polynomial of degree 2. Substituting the 
Legendre expansion for F u and using the orthonormality of the Legendre polynomials and the 
recursion relation (I + l)i^ +1 (/i) = (21 + l)fiPi(fi) — ZP/_i(//), we obtain 

Synchronous gauge — 



Fui 



4 

— ( 
3 

k 2 ( -5, 



2, 



h, 



(7,. 



2(7, 



8 „ 3 r r, 4 i 8 ■ 

15 5 15 5 



2/ + 1 



1) 



Z > 3. 



(49) 



Conformal Newtonian gauge 



1 „ 



(7, 



2/ + 1 



"(l-i) 



/ > 2. 



(50) 



This set of equations governs the evolution of the phase space distribution of massless neutrinos. 
Note that a given mode F\ is coupled only to the (I — 1) and (I + 1) neighboring modes. 



The Boltzmann equation ( |48| ) has been transformed into an infinite hierarchy of moment 
equations that must be truncated at some maximum multipole order /max* 

One simple but 

inaccurate method is to set F v i = for I > / max - The problem with this scheme is that the 
coupling of multipoles in equations and (pO]) leads to the propagation of errors from / max to 
smaller I. Indeed, these errors can propagate to I = in a time r ~ l max /k and then reflect back 
to increasing I, leading to amplification of errors in the "square well" < I < Z max - 

An improved truncation scheme is based on extrapolating the behavior of F u [ to / = l max + 1. 
For I > 1, the multipole moments are gauge-invariant and numerical solutions show that they 
exhibit damped oscillations reminiscent of spherical Bessel functions. In fact, if d T ((p + ip) = 
(choosing conformal Newtonian gauge for simplicity), the time dependence of the exact solution 
of the Boltzmann hierarchy is F v i(k,r) oc ji(kr) for / > 0. If we assume that this relation holds 
approximately for time- varying potentials, using a recurrence relation for spherical Bessel functions 
we get 

(2Zmax "I - 1) 



F 



V (Zmax + 1) 



kr 



F 



win 



F, 



V (l n 



-1) • 



(51) 

we have found 
= 0. However, 



By numerically integrating the Boltzmann hierarchy with different choices of l ma , x , 
that this truncation scheme greatly improves on the simple truncation i^rt max +i) : 
time- variations of the potentials during the radiation-dominated era make even equation ( |5l"| ) a 
poor approximation if Z max is chosen too small. We shall discuss in §8 our choices for Z max . 
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5.4. Massive Neutrinos 

Massive neutrinos also obey the collisionless Boltzmann equation. The evolution of the 
distribution function for massive neutrinos is, however, complicated by their nonzero mass. Prom 
equations (^), the unperturbed energy density and pressure for massive neutrinos (labeled by 
subscripts u h n for HDM) are given by 

Ph = a~ A J q 2 dqdnef (q), P h = V 4 J q 2 dq dSl ^ f (q) , (52) 



where e = e(q, r) = \J q 2 + m 2 a 2 , while the perturbations are 

5p h = a- 4 Jq 2 dqdneMq)*, SP h = V 4 f q 2 dq dQ ^ f (q)^> , 

5T^= a- 4 fq 2 dqdnq ni fo(q)*, Y? hj = a" 4 J q 2 dq dQ ^(mn, - itfy) / (g)¥ . (53) 

Since the comoving energy-momentum relation e(q,r) depends on both the momentum and time, 
we can not simplify the calculations by integrating out the g-dependence in the distribution 
function as we did for the massless neutrinos above (see eq. [46]). Instead of applying equation 



, we expand the perturbation ^ directly in a Legendre series 

oo 

9(k, n, q, t) = ]T(-i) Z (2Z + q, r)P^ ■ ft) . (54) 

1=0 

Then the perturbed energy density, pressure, energy flux, and shear stress in fc-space are given by 

5p h = 4vra" 4 J q 2 dq ef (q)^ , 

SPh = fa- 4 Jq 2 dq^-f (q)* , 
(Ph + Ph)O h = ^ka~ A J q 2 dqqf (q)*i , 

(Ph + P h >h = Y«" 4 /rty/o(g)^2. (55) 

Following the same procedure used for the massless neutrinos, the Boltzmann equation 
becomes 



Synchronous gauge 



= *i + «k-Ji ' 

e b amq 

qk , . ( 1 ■ 2 \ din /o 

= ; ; [^f-i ~ (/ + l)^/+l] , ^>3. 

(2Z + l)e L v ; + J 
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Conformal Newtonian gauge 



^0 = 


qk 
e 


ding 


*1 = 




2* 2 ) - - 


*i = 


qk 




(21 + l)e 



^-V'-ri , (57) 

6q amq 

(Z + l)*,+i], Z>2. 

Because of the ^-dependence in these equations, it requires much more computing time to 
carry out the time integration for the massive neutrino. Bond & Szalay (1983) used a 16-point 
Gauss-Legendre method to approximate the q-integration. We do not use this method and instead 
perform the integration using sixth- or eighth-order Newton-Cotes quadrature (plus a remainder 
obtained by asymptotic expansion) with a g-grid of 128 g-points for every wavenumber k. We 
verified that this was enough to ensure a relative accuracy no worse than 10 -4 by trying the 
integration with twice as many points. Then the perturbations 5ph, 5Ph, &hi an d o~h that enter 



the right-hand side of the Einstein equations are calculated from equations (55) by numerically 
integrating ^o^i and ^ over q. 

As with massless neutrinos, we must truncate the Boltzmann hierarchy for massive neutrinos. 
We have found the following scheme to work well: 

(2W + l)e , . 

Because the higher multipole moments decay rapidly once the neutrinos become nonrelativistic, it 
is possible to choose a much smaller Z max for massive neutrinos than for massless ones. 



5.5. Photons 



Photons evolve differently before and after recombination. Before recombination, photons 
and baryons are tightly coupled, interacting mainly via Thomson scattering (and the electrostatic 
coupling of electrons and ions). In Thomson scatterings, the photon energy hu is assumed to 
be much less than the electron rest mass m e ~ 0.511 MeV and the recoil of the electron in 
the initial electron rest frame is neglected. (We are concerned with the period after neutrino 
decoupling, when T < m e .) The classical differential cross section for Thomson scattering is given 
by da/dQ, = 3or(l + cos 2 #)/16-7r, where ax = 0.6652 x 10~ 24 cm 2 and 9 is the scattering angle 
(e.g., Jackson 1975). After recombination, the universe gradually becomes transparent to radiation 
and photons travel almost freely, although Thomson scattering continues to transfer energy and 
momentum between the photons and the matter. 
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The evolution of the photon distribution function can be treated in a similar way as the 
massless neutrinos, with the exception that the collisional terms on the right-hand side of the 
Boltzmann equation are now present and they depend on polarization. Photons propagating in 
direction n are linearly polarized in the plane perpendicular to h due to scattering of electron 
density perturbations with wavevector k. We shall track both the sum (total intensity) and 
difference (Stokes parameter Q) of the phase space densities in the two polarization states for 
each k and n. We shall denote the former (the momentum-averaged total phase space density 
perturbation, summed over polarizations) by Fj(k,n,r), defined as in equation (f46[). Similarly, 
G~(k,n,r) is the difference of the two linear polarization components. The linearized collision 
operators for Thomson scattering are (Bond & Efstathiou 1984, 1987; Kosowsky 1995) 



an e aT 



c 



1 



F 7 + F l0 + An ■ v e - - (F 7 2 + G 7 o + G 72 ) P2 



fdG, 



an e <JT 



-G-y + - (F 7 2 



+ G70 + G72) (1 — P2 



(59) 
(60) 



V dr / c 

where n e and v e are the proper mean density and velocity of the electrons. The terms proportional 
to P2 come from the polarization-dependence of the Thomson cross section which, when averaged 
over incident directions, give an angular dependence 1 + cos 2 9 in the Thomson cross section even 
for unpolarized radiation. The anisotropic scattering and net polarization were both neglected by 
Peebles & Yu (1970); the former (but not the latter) were included by Wilson & Silk (1980, 1981) 
and most later workers. 



Expanding F^(k,h,r) and Gy(k,h, r) in Legendre series as in equation (46) and using the 



relations h ■ v e 
rewritten as 



-(i6b/k)Pi(k ■ n), F 7 i = 4# 7 /(3fc), and F 7 2 = 2cr 7 , the collision operators can be 



SF 7 



an e <TT 



c 



fdG, 



V dr 



4i 



an e <jT 



c 



( 1 1 \ 

9 h )P x + (&t 7 - -G 70 - -G 72 J P 2 - B-«)'(2I + 1)F 7 ^ 

-. 00 
- (F 7 2 + G 70 + G 72 ) (1 - P2) - EH)'( 2/ + WyiPi 



l>0 



(61) 



(62) 



The left-hand-side of the Boltzmann equation for F 7 and G 7 remain the same as for the massless 
neutrinos, so we obtain 

Synchronous gauge — 

4 



^72 



2o\ 



2 






h, 


~ 3 




5 j — 








15 


? 7 



'1) 5 



3 , n 4 • 8 . 9 1 . 

-fcF 7 3 + — h + -rj - -an e a T o- y + —an e a T (G 7 o + G 72 ) , 
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k 



7 l 



G 



7 l 



21 + 1 

k 

21 + 1 



lF 1 a_ 1) -(l + l)F, 



7(1+1) 



an e a T F^i, Z>3, 



^7(2-1) _ + 1 )C 7 (i+l) 



+ an e o~T 



G~i + - (F 7 2 + G 7 o + G 7 2) (S[q + 



(63) 



Conformal Newtonian gauge 



"3^7 + 4 ^ 



1 



<L — o" 7 I + V + o,n e aT( 



"1) ) 



F 7/ 



S 



9 



-.2 = 2(j 7 = — 7 - -kF j3 - -an e a T a 7 + — an e cr T (G 7 o + G 72 ) , 



k 



G 



7 1 



21 + 1 

k 

21 + 1 



lF 1 a_ 1) -(l + l)F, 



7(1+1) 



an e <JTF>y i , Z > 3 



^ 7 (i_i) - (Z+ l)G 7 (| + i) 



+ an e OT 



1 



G~fi + - {Fj 2 + G y o + G 



7 2J OZO 



r , J 12 

c>/o + -=- 



(64) 



The subscripts 7 and 6 refer to photons and baryons respectively. 



We truncate the photon Boltzmann equations in a manner similar to massless neutrinos 
(eq. |5l|), except that Thomson opacity terms must be added. For I = Z max we replace equations 
I and (|J) by 

1 + L 



F^i 



kF 



7(1-1) 



kG 



T 

l + l 



71 



an e o~TF, 



7i 1 



7(1-1) 



Cw — an e o~TG 



(65) 



5.6. Baryons 

The baryons (and electrons) behave like a non-relativistic fluid described, in the absence 
of coupling to radiation, by the energy-momentum conservation equations ( |29"|) and (|30|) with 
5Pb/5pb = c 2 s = u; <C 1 and c = 0. Since the baryons are very nonrelativistic after neutrino 
decoupling (the period of interest), we may neglect w and 5P/5p in all terms except the acoustic 
term c 2 k 2 5 (which is important for sufficiently high k; note that the shear stress term k 2 a is far 
smaller so we neglect it). Before recombination, however, the coupling of the baryons and the 
photons causes a transfer of momentum and energy between the two components. 

From equation ( |22| ) the momentum density T°j for a given species is related to 8 by 
ik^5T°j = (p + P)9. The momentum transfer into the photon component is represented by 
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an e ax{Ob — 8j) of equations (p3| ) and (|6J). Momentum conservation in Thomson scattering then 
implies that a term (4p 7 /3p&) an e (JT(0j — 0&) has to be added to the equation for b (where we 
have used P b <C pf,), so equations ( p9| ) and ([?(]) are modified to become 

Synchronous gauge — 

4 = ~0b- -^h, 

Ob = --0 6 + c^ 2 5 b + ^an e a T (# 7 -0 fe ), (66) 
a 3p& 



Conformal Newtonian gauge 



h = -Ob + 3<j> , 

b = --Ob + cl^db + ^aneO-T^-Obj + k 2 ^. (67) 
a 3p b 



The square of the baryon sound speed is evaluated from 

2 P b k B T b ( 1 dhiT b 
c a = — = 1 - - 



p b p \ 3 din a 

where p is the mean molecular weight (including free electrons and all ions of H and He) and, 
in the second equality, we have neglected the slow time variation of p. (This approximation is 
adequate because even during recombination, when p is largest, the baryons contribute very little 
to the pressure of the photon-baryon fluid.) The baryon temperature evolves according to 

f b = -2-T b + % — ^an e a T (T 7 - T b ) . (69) 
a 6 m e p b 

We assume that electron-ion collisions are rapid enough for kinetic equilibrium to hold with 
a common temperature T b for electrons and all baryon species. Equation (^) follows from 
the first law of thermodynamics, dQ = (3/2)d(Pi,/ p^) + Pbd(l/p b ), with specific heating rate 
Q = 4(p 7 / p b )an e a T k B {T^ - T b ). 



5.7. Tight-Coupling Approximation 



Before recombination the Thomson opacity is so large that photons and baryons are tightly 
coupled, with an e ax = t~ 3> a/a ~ t . The large values of the Thomson drag terms in 
equations (|63|)-(|67|) for # 7 and 6 b make them numerically difficult to solve. Therefore, in this limit 
we shall follow the method of Peebles & Yu (1970) to obtain an alternative form of the equations 
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valid for t c /t -C 1 and kr c <C 1. The starting point is the exact equation obtained by combining 
the second of equations and (|67|), 



(1 + R)0 b + -0 b - c 2 s k 2 5 b 
a 



k 2 R 



4*r 



ct 7 + #(0 7 - e b ) = (i + R)k z i) , 



(70) 



where the right-hand side is included only in the conformal Newtonian gauge; in the synchronous 
gauge it is set to zero. We have defined R = (4/3)p 7 /p b . We shall see that the terms proportional 
to (# 7 — b ) and <r 7 may be neglected to lowest order in max [kr c , t c /t], with the result that the 
baryons and photons behave like a single coupled fluid with sound speed c pb . (This should be 
distinguished from c s , which denotes the sound speed for the baryons only and not that for the 
coupled photon-baryon fluid.) However, we require a more accurate approximation to account for 
the slip between the photon and baryon fluids. 



From the second of equations (p4|), we have 



Ob — 0~ 



k 



(7- 



(71) 



in the conformal Newtonian gauge; in the synchronous gauge we simply set '0 = 0. Writing (9 7 
Ob + (#7 — Ob) and using equation fl70D, we get 



1 + R 



-0 b + k 2 
a 



1 

4^ 



c 1^b 7^-7 + 



+ 



a result that is valid in both gauges. From the third of equations (B3), we have 



9 V3 



7 3 



as 



(72) 



(73) 



in the synchronous gauge; in the conformal Newtonian gauge one sets h = f] = 0. We see that 



(Try ~ 



5 7 x max [kr c , t c /t] (the case kr c corresponding to acoustic oscillations with # 7 



Higher moments of the photon distribution are smaller by additional powers of kr c and we shall 
neglect them in the limit of tight coupling considered here. Our goal is to obtain equations for Ob 
and 0y that are accurate to second order in r c . 



To get an equation for Ob we differentiate equation (ff2j) and use equation (|70|). Assuming that 
the gas is nearly fully ionized so that n e oc a~ 3 and that the baryon temperature is approximately 
the radiation temperature implying c 2 oc a -1 , we obtain 



2R a 



1 + Ra 



+ 



l + R 



a 
-( 

a 



-k 



1 



(L + + k 2 [ c 2 s 5, 



If 
4^ 



+ 0(r 2 ). (74) 



This equation holds in the conformal Newtonian gauge; in the synchronous gauge one should set 
ip = 0. Note that 5b and 5 7 are to be evaluated using the first of equations (|63|)-(|67|). Substituting 
equation ( |74| ) into equation ( |70| ) yields our desired equation of motion for max [kr c , t c /t] <C 1. If 
this condition is violated, then one should use the explicit form of equations (^) and (37) for Ob- 
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To obtain an equation for # 7 we combine the explicit equations for # 7 and 9 b to obtain the 
exact equation 

<V „2t„2 X \ , 1.2 f 1 , „ \ , + 



# 7 = -fi- 1 + -e b - cjk 2 5 b j + fe \-6 7 - a y j + v R ' k 1 ^ (75) 
in conformal Newtonian gauge; in synchronous gauge one sets ip = 0. For 9 b we use the 



tight-coupling approximation (substituting eq. [74] into eq. ||70||) at early times and the exact 



explicit equations (66) or (|67|) at late times. In practice, we switch to the explicit scheme for b 
when T b = 2 x 10 4 K; we switch to the explicit scheme for F y i for I > 1 and T b = 2 x 10 5 K (at 
earlier times these moments are set to zero). We have verified that these switches occur early 
enough to preserve good accuracy in the resulting photon phase space distribution. 



5.8. Recombination 



In order to compute the Thomson scattering terms in the equations of motion for photons and 
baryons we need to know the free electron density n e {r). Our treatment is based on Peebles (1968; 
see also Peebles 1993) with the addition of helium with mass fraction Y = 0.23. We summarize 
the procedure here. 

At high temperatures, hydrogen and helium are both fully ionized. Because of its much 
larger ionization potentials, helium recombines while hydrogen is still fully ionized and the 
free electron density is substantial. Consequently, the helium recombination rates are much 
larger than the expansion rate until helium recombination is essentially complete, so that Saha 
ionization equilibrium is an excellent approximation. We define the helium ionization fractions 
x\ = n(He + )/n(He) and X2 = n(He ++ )/ra(He), where n(He) is the total number density of helium 
nuclei. The Saha equation is 

n e x n +i _ ^9n+i ( m e k B T b \ 3/2 - X n/k B T b 



g n V 2-nti 2 



e -Xn/*B't f ( 76 ) 



where n = or 1 (xq = 1 — x\ — X2), the statistical weights for helium are go = g\ = 1 and 52 = 2, 
and Xn is the ionization potential from n-times ionized helium. 

When hydrogen recombines, the rapidly declining free electron density leads to a breakdown 
of ionization equilibrium. One must integrate the appropriate kinetic equations. Because helium 
recombination is completed before hydrogen begins to recombine appreciably, it is sufficient 
now to treat the helium as being fully neutral. We define the ionization fraction of hydrogen 
as xh = n e /riH where nn is the total number density of hydrogen nuclei. The ionization rate 
equation is (Peebles 1968; Spitzer 1978) 

dXH 



dr 



aC r 



P(T b )(l-x H )-n H a^(T b )x z H . (77) 
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The factor C r is discussed below. The collisional ionization rate from the ground state is 

m P k B T b \ 3/2 ~B 1 /k^T hn {2) i 



m) = e -*l*M a w {Tb) (78) 

where B\ = m e e 4 /(2?i 2 ) = 13.6 eV is the ground state binding energy, and the recombination rate 
to excited states is 

This expression for ^2(^5) is a good approximation (better than one percent for T b < 6000 K). At 
high temperatures this expression underestimates fa but the neutral fraction is negligible so that 
we make no significant error by setting 2 = for T b > B\/k-Q = 1.58 x 10 5 K. 

Recombination directly to the ground state is inhibited by the large Lyman alpha and Lyman 
continuum opacities. Net recombination must occur either by 2-photon decay from the 2s level, 
with a rate A2 S ^i s = 8.227 s _1 , or by the cosmological redshifting of Lyman alpha photons away 
from the line center. Peebles (1968) gives a detailed discussion of these atomic processes. The net 
recombination rate to the ground state is reduced by the fact that an atom in the n = 2 level 
may be ionized before it decays to the ground state. The reduction factor C r is just the ratio of 
the net decay rate (including 2-photon decay and Lyman alpha production at the rate allowed by 
redshifting of photons out of the line) to the sum of the decay and ionization rates from the n = 2 
level: 

Cr ~ A a + A 2 _ ls + /3(2)(T fe ) ' (8Uj 

where 

pM(T b )=P(T b )e +hUa/k * Tb , A-a = ~ 2"TJT~ — ' A Q = ^ = 1.216 x 10- 5 cm , (81) 

where u a = c/X a . For T b <C 10 5 K, it is a very good approximation to replace the number density 
nis of hydrogen atoms in the Is state by (1 — xh)tih- 

We integrate equation ( |77|) using a stable and accurate semi-implicit method with a 
large number of timesteps through recombination. Since the results are independent of the 
perturbations, we pre-compute the ionization history of a model and later use cubic splines 
interpolation to obtain n e (r) (including electrons from both hydrogen and helium) accurately 
during integration of the perturbation equations. The baryon temperature and sound speed are 
also pre-computed for cubic splines interpolation. 



6. Microwave Background Anisotropy 
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Perturbations in the photon phase space distribution correspond to fluctuations in the cosmic 
microwave background radiation. In this section we show how to compute the anisotropy and its 
power spectrum from the results of integration of the Boltzmann equation. 

The photon brightness temperature perturbation A = AT/T is defined by 

f( x \ q ,n J ,r)=f ^ T ^j , (82) 

where fo(q) is the Bose-Einstein distribution of equation (j33|) with e = q. In general the anisotropy 
is a function of all phase-space variables: A = A(x l ,q,rij,T). From equation (^) we obtain, to 
first order in ^f, 

A= _(^Ay\. (83) 



din q 

Because the q-dependence of both the gravitational source terms and the linearized collision 
operator for Thomson scattering in the Boltzmann equation for ^ (eqs. |^0| and [^]) are both 
proportional to dlnfo/dlnq, A is independent of q. In other words, the perturbed microwave 
background radiation has a Planck spectrum with the brightness temperature at a given position 
depending only on the photon direction. (Scattering by a nonlinearly perturbed medium can, 
however, change the spectrum. An example is the Sunyaev-Zel'dovich effect occurring when the 
electron gas is much hotter than the radiation.) In the absence of spectral distortions, A is related 
very simply to the momentum-averaged phase space density perturbation (i.e., the relative energy 
density perturbation): A = ±F 7 . Note that Efstathiou & Bond (1984, 1987) and Bond (1995) 
define A to be the photon density perturbation, 4 times our definition. 

To compute the anisotropy at a given spacetime point (x*,r) we superpose plane- wave 
contributions: 

A{x,h,T)= f d 3 ke il5; A{k,n,T) = fd 3 ke il£ f^(-i) l (2l + l)Ai(k,T)P l (k-n), (84) 
J 1=0 

where A/ = jFji. The anisotropy at the origin may be expanded in spherical harmonics in the 
usual manner: 

OO Z „ 

A (™) = E E a imYi m (n) , a lm = (-i) 1 4vr / d 3 k Y^(k) A t (k, r) . (85) 

1=0 m=-l 

The angular power spectrum of the anisotropy is defined by the covariance matrix of these 
expansion coefficients: 

(aimal'm') = Q 8ui 5mm' , (86) 

where the angle brackets denote a theoretical ensemble average. The angular power spectrum is 
related to the angular two-point correlation function by 

C(6) = (A(ni) A(n 2 )> = — £(2Z + 1) Q P,(ni ■ n 2 ) (87) 

1=0 
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where cos0 = h\- hi. Note that for I > 1 the anisotropy coefficients A/ are gauge- invariant. 
Excluding the monopole and dipole anisotropy, one obtains the same result for the angular power 
spectrum and angular correlation function in both synchronous and conformal Newtonian gauges. 

The anisotropy coefficients A;(/c,r) are random variables with amplitudes and phases 
depending on the initial perturbations. We assume that the initial conditions can be specified in 
terms of the conformal Newtonian gauge potential ip as described in §7. Because the evolution 
equations for F 7 j(fe,r) (and therefore A;) are independent of k, we may write 

A,(fc,r) =Mk)Mk,r) , (88) 

where ipi(k ) is the initial perturbation and Aj(fc, r) is the solution of the Boltzmann equation with 
ipi = 1 (i.e., it is the photon transfer function). Using the two-point correlation function of -0i in 
Fourier space, 

(Mh) Mh)) = Pf(k) Soih + k 2 ), (89) 
where Sd is the Dirac delta function, we obtain the desired result 

Ci = Att J d 3 kP^(k)Af(k,T) . (90) 

As an illustration of this formalism, consider the anisotropy on large angular scales arising 
from the Sachs- Wolfe effect. For isentropic perturbations on scales larger than the acoustic 
horizon at recombination, the present (r = To) anisotropy is related to the perturbation in ip at 
recombination by A(n, To) ~ \ip{x = —nx,T ICC ) with X = T o — Tree ~ tq (Sachs & Wolfe 1967). 
In this case the radiation transfer function is simply Ai(k,r) = \ji{k\)- Now suppose that the 
primeval spectrum of curvature fluctuations is a power-law, P 1 p(k) = Ax 3 {kx) n ~ 4 oc k n ~ 4 (with 
n = 1 corresponding to equal power on all scales, or the scale-invariant Harrison-Zel'dovich 
spectrum). The resulting angular power spectrum is 



2 n 7r 3 r(3-n)r(^±3pi 



9 F 2 (4^ r (2ZjW| 

For n = 1 this reduces to Q w (8vr 2 /9)A/[l(l + 1)]. 



7. Super-Horizon-Sized Perturbations and Initial Conditions 



The evolution equations derived in the previous sections can be solved numerically once the 
initial perturbations are specified. We start the integration at early times when a given fc-mode 
is still outside the horizon, i.e., kr <C 1 where kr is dimensionless. (We follow common usage in 
referring to waves with kr < 1 as being "outside the horizon" even though r is more appropriately 
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called the comoving Hubble distance.) The behavior of the density fluctuations on scales larger 
than the horizon is gauge-dependent. The fluctuations can appear as growing modes in one 
coordinate system and as constant modes in another. As we will show in this section, this is 
exactly what occurs in the synchronous and the conformal Newtonian gauges. 

We first review the synchronous gauge behavior, which has already been discussed by Press &: 
Vishniac (1980) and Wilson & Silk (1981), although these authors did not include neutrinos. We 
are concerned only with the radiation-dominated era since the numerical integration for all the 
fc-modes of interest will start in this era. At this early time, the massive neutrinos are relativistic, 
and the CDM and the baryons make a negligible contribution to the total energy density of the 
universe: ptotai = Pv + p-y The expansion rate is a/a = t~ 1 . We can analytically extract the 
time-dependence of the metric and density perturbations h, rj, 5, and 6 on super-horizon scales 
(kr <C 1) from equations (|2~l|), ( |49| ) and (|63[). The large Thomson damping terms in equations ( |63| ) 
drive the I > 2 moments of the photon distribution function F~i and the polarization function G~i 
to zero. Similarly, F v \ for I > 3 can be ignored because they are smaller than F U 2 by successive 
powers of kr. Equations ( |21a| ), ( plcj ), (f49|), and ( |63[ ) then give 

T 2 h + rh + 6[(1 - R u )5 y + RJ V ] = , 
4 2 • ■ 1 9 

Su + ^e v + ^h = 0, 8 u -^k\5 u -4a u )=0, (92) 

a v -— (2e v + h + 6fj) = 0, 
15 

where we have defined R u = Pu/ip-y + Pv)- For N u flavors of neutrinos (N v = 3 in the standard 
model), after electron-positron pair annihilation and before the massive neutrinos become 
nonrelativistic, p v /p 1 = (7N v /8) (4/1 1) 4 / 3 is a constant. 

To lowest order in kr, the terms oc k 2 in equations ( |9"2| ) can be dropped, and we have 
9 U = 8^ = 0. Then these equations can be combined into a single fourth-order equation for h: 

d 4 h d 3 h . , 



whose four solutions are power laws: h oc r n with n = 0, 1, 2, and —2. From equations (|9^) we 
also obtain 

h = A + B(kr)- 2 + C{kT) 2 + D(kr), 

S = (1 - R^ + R v 6„ = -lB{kT)~ 2 -\c{kr) 2 - X -D{kr) , 

3 3 fa 

9 = (I - R u )0 y + R U 9 U = —Dk, (94) 

o 

and A, B, C, and D are arbitrary dimensionless constants. The other metric perturbation r\ can 
be found from equation ( |21a| ): 

r, = 2C + ^D(kT)- 1 . (95) 
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Press & Vishniac (1980) derived a general expression for the time dependence of the four 
eigenmodes. They showed that of these four modes, the first two (proportional to A and B) are 
gauge modes that can be eliminated by a suitable coordinate transformation. The latter two 
modes (proportional to C and D) correspond to physical modes of density perturbations on scales 
larger than the Hubble distance in the radiation-dominated era. Both physical modes appear 
as growing modes in the synchronous gauge, but the C(kr) 2 mode dominates at later times. In 
fact, the mode proportional to D in the radiation-dominated era decays in the matter-dominated 
era (Ratra 1988). We choose our initial conditions so that only the fastest-growing physical 
mode is present (this is appropriate for perturbations created in the early universe), in which 
case 7 = 9 U = f\ = to lowest order in kr. To get nonzero starting values we must use the full 



equations (92) to obtain higher order terms for these variables. To get the perturbations in the 
baryons we impose the condition of constant entropy per baryon. Using all of these inputs, we 
obtain the leading-order behavior of super-horizon-sized perturbations in the synchronous gauge: 



Synchronous gauge 



<5 7 = --C(kr) 2 , 5 C = 5 b = j5 u = |j<5 7 , 



e c = o, 0i = 0b = _l_c(k±T 3 ), i5*^ ftr, ( 96 ) 

a - = 3(15 4 + V) (fcT)2 ' 

k = C(krf, V^C-^^Cikrf. 

The initial conditions for the moments I > 1, of the massive neutrino distribution can be related 
to an d the variables above by equations (]56|). To obtain the initial ^o> we can write the perturbed 
distribution function using equation (f§) as / = f (q)(l + *o) = 2/ip 3 {exp[g//c(T + AT)] + l}" 1 , 
where AT/T = 5 U /A by the isentropic condition. Then using equations (|56|) , we find the first three 
moments to be 

1 . d In /o 

6qk a In q 
2 amq 

where the terms proportional to m^a 2 /q 2 in ^ are dropped. The higher moments ^ (I > 3) are 
negligible for kr <C 1. 

The initial conditions for the isentropic perturbations in the conformal Newtonian gauge can 
be obtained either by solving equations (^) , ([5(]) , and (^) , or using the transformations given by 



equations (18) and (p7|) (which enables us to relate the amplitudes in the two gauges). We find for 



the growing mode 
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Conformal Newtonian gauge 



= —2tp , 5 r = 5h = —o~» = — ( 

15 + 4R, ^ 4 4 



7 = ft, = 0C = 9 b = 15 ^ (fe 2 r) = ^V)^ , (98) 
4/7 1 

The massive neutrino moments ^ in this gauge are related to 5 U , 6 U , and a v by the same equations 



(eq. [97]) as in the synchronous gauge. As claimed earlier, tp = <f> to zeroth order in kr when 
no neutrinos are present (i.e., R v = 0). If we characterize the perturbations in the conformal 
Newtonian gauge by the potential tp, all matter and metric variables have a very simple form 
outside the horizon. The neutrino energy fraction R v enters only in the second potential ^ as a 
result of the shear stress produced by the free-streaming neutrinos. Bardeen (1980) was concerned 
that a large shear stress would lead to large metric perturbations in the conformal Newtonian 
gauge. We see that this does not happen for isentropic growing-mode perturbations in which the 
shear stress arises solely due to the free-streaming of relativistic collisionless particles. 

We see that 5 grows with time in the synchronous gauge but remains a constant in the 
conformal Newtonian gauge before horizon crossing. Another significant difference is the larger 
value of the velocity perturbations for small kr in the conformal Newtonian gauge. Physically, this 
difference arises because velocity perturbations vanish to lowest order in the synchronous gauge 
because the synchronous gauge spatial coordinates are Lagrangian coordinates for freely-falling 
observers (§2). The next-order velocity perturbations differ for the neutrinos and photons because 
these two fluids have effectively different equations of state: the neutrinos are collisionless while 
the photons behave like a perfect fluid due to their strong coupling to the baryons. In the 
conformal Newtonian gauge, the lowest-order velocity perturbations do not vanish because the 
conformal Newtonian gauge spatial coordinates are Eulerian coordinates. If we were to include the 
next-order corrections to 8 proportional to fe 4 T 3 , differences between the different fluid components 
would appear in equations (|9§|). 

In the conformal Newtonian gauge, the mode proportional to D in the synchronous gauge 
yields <p tx tp oc 6 oc (fcr)~ 3 . Thus, this mode corresponds to a decaying mode in the conformal 
Newtonian gauge even though it yields S oc (kr) in the synchronous gauge. The two gauge modes 
(A and B in eq. |94|) do not exist in the conformal Newtonian gauge. 



8. Integration Results in Two CDM+HDM Models 
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We apply the results derived in the previous sections to two spatially flat cosmological 
models consisting of a mixture of CDM and HDM: (1) the "old" model with a neutrino fraction 
of Qj, = 0.3 and a corresponding neutrino mass of m v = 93.13 (Q u h 2 ) eV = 7 eV (e.g. Davis et 
al 1992; Klypin et al 1993; Jing et al 1994; Cen & Ostriker 1994), and (2) the tt v = 0.2 model 
{m u = 4.7 eV), which gives a better match to high-redshift observations (Ma & Bertschinger 
1994b; Klypin et al 1995). Both models assume ^baryon = 0.05 and H$ = 50 km s _1 Mpc -1 . 

In Fourier space, all the k modes in the linearized Einstein, Boltzmann, and fluid equations 
evolve independently; thus the equations can be solved for one value of A; at a time. Moreover, 
all modes with the same k (the magnitude of the comoving wavevector) obey the same evolution 
equations. For purposes of computing matter transfer functions, we integrated the equations of 
motion numerically over the range 0.01 Mpc" 1 < k < 10 Mpc" 1 using 31 points evenly spaced 
in log 10 k with an interval of Alog 10 k = 0.1. For computing microwave background anisotropy 
we used linear spacing in k with 5000 points up to k = 0.5 Mpc -1 . The time integration was 
performed using the standard fifth- and sixth-order Runge-Kutta integrator dverk (obtained from 
netlib@ornl.gov). It began at conformal time tq = 3 x 10 -4 Mpc (z ~ 10 9 ), which was chosen 
so that the largest k (i.e. the smallest wavelength) was well outside the horizon at the onset of 
the integration. The full integration was carried to z = 0. The microwave background anisotropy 
calculations were performed on a Cray C-90 at the Pittsburgh Supercomputing Center, using a 
vectorized code that runs at about 570 MFlops per processor. 

The Einstein equations provide redundant equations for the evolution of the metric 
perturbations. In the synchronous gauge we chose to use ah and r] as the primary metric 
perturbation variables in the integration, and used equation ( |21b| ) and a combination of ( |21a| ) 
and ( |21cj ) as the evolution equations. In the conformal Newtonian gauge we integrated (/> using 
equation ( |23b| ) and obtained tjj algebraically using ( 23d] ) . In both gauges we used the time-time 



Einstein equation (eqs. |21a| 1 and [23a]) to check integration accuracy. In the conformal Newtonian 
gauge it is possible to avoid integration of the metric perturbations altogether by combining 
equations fl23a| ) and ( |23b| ) into an algebraic equation for 0. However, we found that this gave 
numerical difficulties because the initial value of 4> has to be set with exquisite precision when 
kr <C 1. We also found it necessary to obtain the initial #'s from <f> and the <5s in the combined 
constraint equations of fl23a|) and ( |23b| ) . Although the analytical expressions in equations (^) are 
good approximations for kr <C 1, slight deviations from the energy- momentum constraints was 
found to cause numerical difficulties. 

In the computations of the potential and the density fields (shown in Figs. [2| - |j), the 
photon and the massless neutrino phase space distributions were expanded in Legendre series 
(see eq. |46|]) with up to 2000 /-values in order to guarantee sufficient angular resolution. In the 
calculations of the temperature fluctuations (Fig. |l]), we used the criterion Z max = 1.5/cr max + 10, 
where a(r max ) = 1 (r max 12000). The massive neutrinos are computationally expensive due to 
the momentum-dependence in equations ( |56| ) and (|57|). For all computations we performed the 
massive neutrino calculations on a grid of 128 (/-points including 50 /-values for every q. Using our 
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truncation schemes given by equations (|51~|), (|58|), and (^), truncating the Boltzmann hierarchies 
at £ max = 50 for massive neutrinos and 2000 for the massless particles was adequate at better than 
the 0.1% level. We checked that all of our numerical approximations are adequate by increasing 
the grids of k, q, and I values as well as decreasing the integration timestep. We estimate that our 
final results have a relative accuracy better than 10 -3 . 

Figure [l| shows the angular power spectrum C\ (defined in §6) of the photon anisotropy in the 
CDM and the CDM+HDM models. We have assumed a scale-invariant spectrum of the primeval 
potential tp normalized to P^(k) = /c~ 3 . Integrations using conformal Newtonian and synchronous 
gauges agree to better than 0.1%. Our results for CDM agree well with those of other workers 
who included accurate numerical integrations with massless neutrinos, polarization, and helium 
recombination (P. Steinhardt 1994, private communication; Hu et al. 1995). To our knowledge, 
ours are the first results for the CDM+HDM models including all of this physics and more as 
described in the preceding sections. 

Inclusion of massive neutrinos increases the anisotropy at large I (by about 10 percent at the 
second and third acoustic peaks; see Fig. |l]) because the decreased perturbation growth leads to 
a (f) contribution to the photon energy density fluctuation growth in equation fl64|). At the same 
time, the smaller <j> at early times with massive neutrinos leads to a slight reduction in C\ below 
the first acoustic peak. The differences between the Q v = 0.2 and 0.3 models are small (middle 
panel of Fig. |l|). Polarization decreases the anisotropy at high I by increasing the photon shear 
stress oVy, leading to increased radiative diffusion (bottom panel of Fig. |l|). 

Figure ^ shows the time evolution of the metric perturbations <j>(k, r) and ip(k, r) in the 
conformal Newtonian gauge for all 31 values of k in the Q u = 0.2 CDM+HDM model. (The 
metric perturbations in the synchronous gauge have no simple physical interpretation, so we shall 
not bother presenting them.) The overall normalization was chosen arbitrarily (corresponding to 
C = — 1/6 in eqs. |M| and |p8[ ). The difference between ip and <p in the radiation-dominated era is 
due to the shear stress contributed by the relativistic neutrinos (eq. |)8| ) which make up a fraction 
R u = 0.4052 of the total energy density. On scales much smaller than the horizon, ip corresponds 
to the Newtonian gravitational potential and <fi = tp in the matter-dominated era. As is well 
known, the potential is constant for growing-mode density perturbations of CDM in an Einstein-de 
Sitter universe. In a mixed dark matter model, however, the CDM density perturbation growth 
can be suppressed by the lack of growth of HDM perturbations so that ip decays slowly. 

The behavior of the metric perturbations can be understood as follows. All of the 31 
fe-modes are outside the horizon at early times when r < 0.01 Mpc. The horizon eventually 
"catches up" and a given /c-mode crosses inside the horizon when kr is about ir. The modes 
with larger k (i.e., shorter wavelengths) enter the horizon earlier. If a given fe-mode enters the 
horizon during the radiation-dominated era, the tight coupling between photons and baryons due 
to Thomson scattering induces damped acoustic oscillations in the conformal Newtonian gauge 
metric perturbations, which are exhibited in Figure [2] by the modes with k > 0.1 Mpc -1 . (In 
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fact, it is not the speed-of-light horizon that sets the scale for the oscillation and damping of the 
potential. Rather, as we show below, it is the acoustic horizon. These horizons are similar during 
the radiation-dominated era because the sound speed of the photon-baryon fluid c p b is c/y/3.) The 
modes with k < 0.1 Mpc -1 enter the horizon during the matter-dominated era and do not oscillate 
acoustically because the Jeans wavenumber kj = {AnxGpa 2 / 'c^) 1 / 2 has then become much larger 
than the wavenumbers under investigation. 

We can understand the damped oscillations more quantitatively by studying the Einstein 
equations ( |23"|) in the conformal Newtonian gauge. Analytical solutions can be found for a perfect 
fluid with no shear stress, in which case (j) = ip. Using c 2 b = 5P/5p ~ p/p and a/a = 2T~ 1 /(l + 3c 2 )fe ) 



(from eqs. [19| and [pQ]), equations ([23;) can be combined to yield 



T 



2 7 , 6(1 + ^ 



+ 1 , q 2 T( t> + (kCpbTTt = , (99) 
1 + 6C pb 



whose solutions (approximating c„b as a constant) are Bessel functions with a power-law pre-factor: 

.-,2 



<j>± = (kc pb T) U J ±u {kc ph T) , V = ,o V 2\ ■ ( 100 ) 

Z \ l + 6c pb> 



The <p- solution corresponds to the decaying mode discussed in §6, so we shall ignore it. In the 

c lb = s aud i/ = §, 



radiation-dominated era, ci b = i and v = |, so that 



0/2 f constant, k(,T«l, 

4> + = {kc ph r) ^3/2(^r)oc| a _ 2cos( ^ r)) kCpbT>>1 (101) 

We refer to the damping for kc p bT > 1 as acoustic damping; it corresponds to constant-amplitude 
acoustic oscillations in the density contrast of the photon-baryon fluid. The damped oscillations 
are apparent in Figure ^ for r < r eq . The analytic solution holds of course only in the absence of 
neutrinos; obtaining the correct amplitudes for ip and (f) in CDM+HDM models shown in Figure || 
required the full integration discussed in this paper. 

After the universe becomes matter-dominated, acoustic damping of the potential ceases, 
and the only physical process causing the potential to change is the free-streaming damping of 
perturbations in the massive neutrinos. We shall discuss this further after examining the evolution 
of the density perturbations. 

Figure || shows the evolution of the density perturbations for the five particle species in 
the £l u = 0.2 CDM+HDM model in the two gauges from our numerical integration. Three 
wavenumbers are plotted: k = 0.01 Mpc^ 1 (Fig. |a), k = 0.1 Mpc^ 1 (Fig. |b), and k = 1.0 Mpc^ 1 
(Fig. [|c). Each mode is normalized with the same initial amplitude for cf> as in Figure ^. There 
are several notable features: 

Before horizon crossing — 

(1) The initial amplitudes of the <5's are related by the isentropic initial conditions: 
<5 7 = 5 U = 5h = 4:5b /3 = 4<5 c /3. The behavior of 5 outside the horizon is strongly gauge-dependent. 
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In the synchronous gauge, Figure [3| shows that all the <5's before horizon crossing in the 
radiation-dominated era grow as a 2 . This confirms equations ( |96| ) since a(r) oc r at this time. It 
is straightforward to show that in the matter-dominated era (for 0=1), 5 oc r 2 oc a for all modes 
before horizon crossing. In the conformal Newtonian gauge, the 5's remain constant outside the 
horizon as derived in equations (|98|). 

After horizon crossing — 

(2) The perturbations come into causal contact after horizon crossing and become nearly 
independent of the coordinate choices. As Figure || shows, 5 C , 6b, and 5h in the two gauges 
are almost identical at late times. For CDM, the /c-modes that enter the horizon during the 
radiation-dominated era behave very differently from those entering in the matter-dominated era. 
The critical scale separating the two is the horizon distance at the epoch of radiation-matter 
equality (a eq ~ 2 x 10~ 4 f2/i 2 ): k cq = 2ir/T eq ~ 0.1 Mpc -1 for our parameters. For the modes 
with k > 0.1, horizon crossing occurs when the energy density of the universe is dominated by 
radiation; thus the fluctuations in the CDM can not grow appreciably during this time. For the 
photons and the baryons, the important scale is the horizon size at recombination (a rec ~ 10~ 3 ): 
k rec = 2it/t tcc ~ 0.025 Mpc -1 . The modes with k > 0.025 (see Figs. |^b and ||c) enter the horizon 
before recombination, so the photons (long-dashed curves) and baryons (dash-dotted curves) 
oscillate acoustically while they are coupled by Thomson scattering. The coupling is not perfect. 
The friction of the photons dragging against the baryons leads to Silk damping (Silk 1968), which 
is prominent in Figure |c at a ~ 10 -3 ' 5 . The baryons decouple from the photons at recombination 
and then fall very quickly into the potential wells formed around the CDM, resulting in the rapid 
growth of 5b in Figures |b and |c. 

(3) Neutrinos decouple from other species at T ~ 1 MeV and a ~ 10~ 10 . At this early time, 
both the massless and the eV-range massive neutrinos behave like relativistic collisionless particles. 
The massive neutrinos become non-relativistic when 3/cbT„ ~ m„, corresponding to a nr ~ 10 -4 
for 4.7 eV neutrinos. Close inspection of the figures at a ~ a nr reveals that 5^ (short-dashed 
curve) is indeed making a gradual transition from the upper line for the radiation fields to the 
lower line for the matter fields. Although the Jeans length of a fluid is not well defined for 
collisionless particles such as the neutrinos, the criterion for free-streaming damping is similar 
to the Jeans criterion for gravitational stability: free-streaming is important for k > kf s , where 
kf s (a) = 47rG / oo 2 /t; 2 ned and v me d is the median neutrino speed. When the neutrinos are relativistic, 
v ~ 1 and ki s (a) oc a -1 for a < a eq . After the neutrinos become non-relativistic, the median 
neutrino speed is v me d = 'ik^T^^j am v = 15a _1 (m iy /10eV) _1 km s . In the matter-dominated 
era, we have AnGpa? = ^HqCi~ 1 from the Friedmann equation, and therefore 

k is (a) = 8a 1 / 2 hMpc- 1 . (102) 

In Figure |3)a, since horizon crossing occurs when the free-streaming effect is already unimportant, 
the evolution of 5h is very similar to that of CDM. In Figures [|b and ||c, however, the free 
streaming effect is evident and the growth of 5h is suppressed until k{ s (a) grows to ~ k. After 
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kf s (a) > k, 5h can grow again and catch up to 5 C . Since kf s oc a 1 / 2 , the larger k modes suffer 
more free-streaming damping and 5h can not grow until later times. The damping in 5h also 
affects the growth of 5 C , slowing it down more for larger Q u compared to the pure CDM model. 
This effect is apparent in the z = power spectra of the pure CDM and the two CDM+HDM 
models in Figure Contrary to previous figures, the curves in Figure || are normalized to the 
COBE-compatible rms quadrupole moment of Qrms-PS = 17.6 (Bennett et al. 1994). 



9. Summary 

The purpose of this paper was to present in both the synchronous and the conformal 
Newtonian gauges a complete discussion of the linear theory of scalar gravitational perturbations 
that is applicable to any flat CDM, HDM or CDM+HDM model (including a possible cosmological 
constant). For historical reasons, most calculations of linear fluctuation growth have been carried 
out in the synchronous gauge. The conformal Newtonian gauge, however, offers an alternative 
that is free of the gauge ambiguities and coordinate singularities associated with the synchronous 
gauge. We derived the coordinate transformation relating the two gauges and presented in parallel 
in both gauges the complete set of evolution equations: the Einstein equations for the metric 
perturbations, the Boltzmann equations for the photon and neutrino phase space distributions, 
and the fluid equations for CDM and baryons. A detailed discussion of the microwave background 
anisotropy calculations was also presented. Care was taken to include all important higher 
moments of the neutrino phase space distribution and the effects of helium recombination and 
photon polarization. 

We solved the linear theory in the standard CDM model and two spatially flat CDM+HDM 
models with Q u = 0.2 and 0.3, assuming a scale-invariant spectrum of isentropic primordial 
fluctuations. (The baryon fraction was taken to be ^baryon = 0-05 and Hq=50 km s _1 Mpc -1 .) 
The evolution of the metric perturbations and the density fields for all five particle species were 
presented, along with the first accurate calculations of the photon anisotropy power spectrum in 
CDM+HDM models. We also illustrated the gauge dependence of the density fields before horizon 
crossing and discussed the physical interpretation of the results. 

Interested users may obtain our programs to integrate the perturbation equations at 
trttp : //arcturus . mit . edu/ cosmics. 
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10. Figure Captions 

Fig. [l]: Top: The angular power spectrum C\ of the photon anisotropy (including polarization) in 
the CDM (solid), the Vt v = 0.2 CDM+HDM (dashed), and the Q u = 0.3 CDM+HDM (dotted) 
models. Middle: The ratio of Q for the CDM+HDM models relative to the CDM (dashed for 
Q. v = 0.2; dotted for VL U = 0.3). Bottom: The fractional difference in C% with and without 
polarization for the three models. 

Fig. [2]: The scalar metric perturbations (j){k } r) and ip(k,r) in the conformal Newtonian gauge in 
the Sly = 0.2 CDM+HDM model. The 31 curves from left to right correspond to 31 values of 
k between 10.0 Mpc -1 and 0.01 Mpc -1 . The labels r nr , r eq and r rec indicate, respectively, the 
time the 4.7 eV neutrinos become non-relativistic, the matter-radiation equality time, and the 
recombination time. 

Fig. ||: Evolution of the density fields in the synchronous gauge (top panels) and the conformal 
Newtonian gauge (bottom panels) in the Vt u = 0.2 CDM+HDM model for 3 wavenumbers k= 
0.01 (Fig. ^a), 0.1 (Fig. |3]b) and 1.0 (Fig. ||c) Mpc -1 . In each figure, the five lines represent 
5 C , 5b, <5 7 , S u and Sh for the CDM (solid), baryon (dash-dotted), photon (long-dashed), massless 
neutrino (dotted), and massive neutrino (short-dashed) components, respectively. 

Fig. f§ The z = power spectra for the pure CDM (solid), the Q v = 0.2 CDM+HDM (dashed), 
and the Q v = 0.3 CDM+HDM (dotted) models. In the mixed models, the upper curve is for the 
CDM component and the lower one is for the HDM. The COBE-compatible quadrupole moment 
of Qrms-PS = 17.6 fj, K is used. 
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Fig. 1. — Top: The angular power spectrum C\ of the photon anisotropy (including polarization) 
in the CDM (solid), the tt u = 0.2 CDM+HDM (dashed), and the Sl u = 0.3 CDM+HDM (dotted) 
models. Middle: The ratio of C/ for the CDM+HDM models relative to the CDM (dashed for 
£l u = 0.2; dotted for £l u = 0.3). Bottom: The fractional difference in C/ with and without 
polarization for the three models. 
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Fig. 2. — The scalar metric perturbations (f)(k, r) and ip(k, r) in the conformal Newtonian gauge 
in the Q. v = 0.2 CDM+HDM model. The 31 curves from left to right correspond to 31 values 
of k between 10.0 Mpc" 1 and 0.01 Mpc" 1 . The labels 

Tnr; T"eq and T re c indicate, respectively, the 
time the 4.7 eV neutrinos become non-relativistic, the matter-radiation equality time, and the 
recombination time. 
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Fig. 3. — Evolution of the density fields in the synchronous gauge (top panels) and the conformal 
Newtonian gauge (bottom panels) in the Cl v = 0.2 CDM+HDM model for 3 wavenumbers k= 
0.01 (Fig. ||a), 0.1 (Fig. ^b) and 1.0 (Fig. ||c) Mpc -1 . In each figure, the five lines represent 
5 C , Sbfd^jSu and 8h for the CDM (solid), baryon (dash-dotted), photon (long-dashed), massless 
neutrino (dotted), and massive neutrino (short-dashed) components, respectively. 
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Fig. 4.— The z = power spectra for the pure CDM (solid), the 9. v = 0.2 CDM+HDM (dashed), 
and the Q. v = 0.3 CDM+HDM (dotted) models. In the mixed models, the upper curve is for the 
CDM component and the lower one is for the HDM. The COBE-compatible quadrupole moment 
of Q 

rms— PS = 17.6 fiK is used. 



